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ABSTRACT 

We develop a sub-grid model for the growth of supermassive Black Holes (BHs) and 
their associated Active Galactic Nuclei (AGN) feedback in hydrodynamical cosmologi- 
cal simulations. This model transposes previous attempts to describe BH accretion and 
AGN feedback with the Smoothed Particle Hydrodynamics (SPH) technique to the 
Adaptive Mesh Refinement (AMR) framework. It also furthers their development by 
implementing a new jet-like outflow treatment of the AGN feedback which we combine 
with the heating mode traditionally used in the SPH approach. Thus our approach 
allows one to test the robustness of the conclusions derived from simulating the impact 
of self-regulated AGN feedback on galaxy formation vis-a-vis the numerical method. 
Assuming that BHs are created in the early stages of galaxy formation, they grow by 
mergers and accretion of gas at a Eddington-limited Bondi accretion rate. However 
this growth is regulated by AGN feedback which we model using two different modes: 
a quasar-heating mode when accretion rates onto the BHs are comparable to the Ed- 
dington rate, and a radio-jet mode at lower accretion rates which not only deposits 
energy, but also mass and momentum on the grid. In other words, our feedback model 
deposits energy as a succession of thermal bursts and jet outflows depending on the 
properties of the gas surrounding the BHs. We assess the plausibility of such a model 
by comparing our results to observational measurements of the coevolution of BHs 
and their host galaxy properties, and check their robustness with respect to numeri- 
cal resolution. We show that AGN feedback must be a crucial physical ingredient for 
the formation of massive galaxies as it appears to be able to efficiently prevent the 
accumulation of and/or expel cold gas out of halos/galaxies and significantly suppress 
star formation. Our model predicts that the relationship between BHs and their host 
galaxy mass evolves as a function of redshift, because of the vigorous accretion of 
cold material in the early Universe that drives Eddington-limited accretion onto BHs. 
Quasar activity is also enhanced at high redshift. However, as structures grow in mass 
and lose their cold material through star formation and efficient BH feedback ejection, 
the AGN activity in the low-redshift Universe becomes more and more dominated by 
the radio mode, which powers jets through the hot circum-galactic medium. 

Key words: galaxies: evolution - galaxies: active - galaxies: jets - galaxies: quasars: 
general - methods: numerical 
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1 INTRODUCTION 

Evidence for the ubiquitous presence of supermassive black 
holes (BH) in the centres of galaxies is overwhelming 
(Kormendy & Richstone 1995). BHs spanning a range of 
masses from a few 10^ M0 in the centre of galaxies with 
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small bulges like our Miky-Way (Schodel 2002) up to sev- 
eral 10® Mq for elliptical galaxies in the cores of groups and 
clusters of galaxies (Magorrian et al. 1998) have now been 
observed. These supermassive BHs are not only seen in the 
near universe, but very luminous quasars discovered beyond 
z > 6 (Fan 2003) suggest that they are already in place dur- 
ing the early stages of galaxy formation. As a consequence, 
it is now widely accepted that a large variety of galaxies host 
BHs in their centres, and that these BHs somehow influence 
the evolution of their host galaxies. 

Observations by Magorrian et al. (1998) first pointed 
out a relationship between the central BHs and their 
host galaxy bulge mass with a quasi-linear scaling 
(Laor 2001; McLure & Dunlop 2002; Marconi & Hunt 2003; 
Haring & Rix 2004). A similar, albeit arguably tighter 
correlation is also found between the BH mass and 
the stellar velocity dispersion (Ferrarese & Merritt 2000; 
Gebhardt et al. 2000; Tremaine et al. 2002; Giiltekin et al. 
2009), or the Sersic index that measures the concentra- 
tion of the bulge (Graham & Driver 2007). These corre- 
lations define a BH fundamental plane, similar to the 
fundamental plane of elliptical galaxies, that links BHs 
to bulge stellar masses, velocity dispersions and effective 
radii (Hopkins et al. 2007). 

These observations led to the suggestion that the 
growth of BHs is self-regulated by the energy released 
during their accretion phase. This would be sufficient to 
unbind the gas in the galaxy and form powerful out- 
flows (Silk & Rees 1998; King 2003; Wyithe & Loeb 2003). 
There exist abundant observational evidence of such out- 
flows with direct imaging of X-ray cavities in the vicinity of 
elliptical galaxies (Boehringer et al. 1993; Owen et al. 2000; 
Birzan et al. 2004; McNamara et al. 2005; Fabian et al. 
2006; Taylor et al. 2006; Dong et al. 2010; Dunn et al. 2010) 
or indirect detection of Broad Line Absorption regions in 
the spectra of quasars (Chartas et al. 2003; Crenshaw et al. 
2003; Pounds et al. 2003). These observations are sup- 
ported by numerical models of the microphysics of BH ac- 
cretion discs that can drive massive hydro-magnetic out- 
flows and jets (De Villiers et al. 2005; McKinney 2006; 
McKinney & Blandford 2009), and large amounts of heat 
carried by photons that could potentially ionize the sur- 
rounding gas. Feedback from BHs is commonly called Active 
Galactic Nuclei (AGN) feedback as the energy is emitted 
from the centres of galaxies where BH reside. However, this 
generic appellation encompasses various modes of energy re- 
lease from the central source. 

It is commonly believed that two of these modes can 
describe AGN feedback. These are very similar to the two 
modes observed for X-ray binaries (Churazov et al. 2005; 
Merloni & Heinz 2008). The so-called 'radio' mode is asso- 
ciated with a strong radio emission filling X-ray depressed 
cavities with relativistic electrons. Most of the energy in the 
radio mode is driven by a mechanical jet feedback mech- 
anism infiating the radio lobe itself. This mode is equiv- 
alent to the 'low/hard' state of X-ray binaries that has a 
hard X-ray power spectrum with a cut-off at a few 100 
keV, and gas accretion at low Eddington rates. The same 
tendency has been confirmed for supermassive BHs, for 
which the radio loudness is stronger at lower Eddington 
accretion ratios (Nagar et al. 2005; Chiaberge et al. 2005; 
Churazov et al. 2005). This mode is clearly associated with 



mechanical feedback where most of the energy is powered by 
the jet mechanism and largely overwhelms the X-ray con- 
tribution from the nucleus, as in the well-studied case of 
M87 (Owen et al. 2000). 

A transition to a radio-quiet 'quasar' mode occurs at a 
few ~ 10^^ of the Eddington accretion rate for X-ray bina- 
ries (Maccarone 2003). Above this threshold. X-ray binaries 
enter a 'high/soft' state emitting a soft and thermal X-ray 
spectrum with almost no trace of a jet mechanism, which is 
the equivalent of a quasar spectrum. The thermal emission is 
well described by the standard model of optically thick and 
geometrically thin accretion discs from Shakura & Sunyaev 
(1973), whereas the launching of the jet for the ra- 
dio mode comes from optically thin and geometrically 
thick (and radiatively inefficient) accretion discs modeled 
with the Advection-Dominated Accretion Flow (ADAF) 
from Narayan & Yi (1994) or with the Adiabatic Inflow- 
Outfiow Solutions (ADIOS) from Blandford & Begelman 
(1999). 

On the other hand, observations of tidally dis- 
rupted galaxies often reveal powerful AGN activ- 
ity (Barnes & Hernquist 1992). Mergers between galaxies 
are invoked to compress the inter-stellar medium (ISM) 
and provide a fresh flow of material onto the BH in that 
case. Theoretical models of galaxy formation that associate 
the growth of BHs to such events have been successful at 
reproducing many properties of the population of quasars as 
well as the BH density seen at low redshift (Cattaneo et al. 
1999; Kauffmann & Haehnelt 2000; Granato et al. 2001; 
Volonteri et al. 2003). Semi-analytic models of galaxy 
formation require AGN feedback to suppress the cooling 
catastrophe in massive galaxies, match the bright-end 
of the galaxy luminosity function, and obtain bulge- 
dominated galaxies (Croton et al. 2006; Cattaneo et al. 
2006; Bower et al. 2006, 2008; Somerville et al. 2008). 

Important steps have recently been incorporated with 
hydrodynamical simulation to better quantify the neg- 
ative AGN feedback effect on star formation and gas 
properties of galaxies. These involved simulating ideal- 
ized disc galaxies and galaxy mergers (Springel et al. 2005; 
Nayakshin & Power 2010; Debuhr et al. 2010), or clus- 
ters of galaxies (Cattaneo & Teyssier 2007; Dubois et al. 
2009). Self-consistent sub-grid models of the AGN feed- 
back heating mode have been introduced in ACDM cos- 
mological simulations using the Smoothed Particle Hydro- 
dynamics (SPH) technique as implemented in the GADGET 
code (Sijacki et al. 2007; Booth & Schaye 2009). Alterna- 
tive approaches based on the injection of Cosmic Rays have 
also been explored (Sijacki et al. 2008). These cosmologi- 
cal simulations have successfully reproduced relationships 
between BH masses and galaxy properties (Sijacki et al. 
2007; Di Matteo et al. 2008; Booth & Schaye 2009, 2011), 
and suppressed the cooling catastrophe in groups and clus- 
ters of galaxies (Puchwein et al. 2008; Khalatyan et al. 2008; 
McCarthy et al. 2010, 2011). Recently, AGN feedback asso- 
ciated to BH growth has been introduced in Adaptive Mesh 
Refinement (AMR) cosmological re-simulations of a galaxy 
cluster with the RAMSES code. These featured either a jet- 
kinetic mode (Dubois et al. 2010) or a thermal energy in- 
put (Teyssier et al. 2011) (see also Dubois et al. 2011), but 
the cosmic co-evolution of BHs and galaxies has not yet been 
studied using grid based techniques. 
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As limitations of the standard SPH technique to cap- 
ture Kelvin-Helmholtz instabihties have been pointed out 
(Agertz et al. 2007), and authors like Mitchell et al. (2009) 
have shown that it could have severe consequences on the 
properties of the intra-cluster gas where AGN-host ellipti- 
cals reside, it seems worthwhile to investigate this important 
issue using a different numerical technique. Therefore, this 
work uses an Eulerian grid based approach to model AGN 
feedback from BHs by a dual jet/heating sub-grid model 
representative of the radio and quasar mechanisms with the 
self-regulated growth of BHs. We emphasize the importance 
of these two different modes on the long-term evolution of 
galaxies, trying to explain why quasars are a common in- 
gredient of the young Universe and why radio jets are more 
commonly observed in late and massive structures such as 
clusters of galaxies. 

The paper is organised as followed. In section 2, we 
detail the numerical technique used for following the BH 
growth along with its associated AGN feedback, and the 
standard models employed for galaxy formation (cooling, 
star formation, supernovae feedback, etc.). In section 3, we 
describe the set of simulations employed to test the AGN 
feedback model with the RAMSES code. In section 4, we 
present a parameter study of the AGN feedback model and 
Eissess the convergence of the results vis-a-vis resolution. Sec- 
tion 5 scrutinizes what drives the domination of the quasar 
and the radio mode at different epochs. Finally, we comment 
on the results in section 6. 



2 MODELLING THE PHYSICS OF GALAXY 
FORMATION 

2.1 BH growth and AGN feedback 

Sink particles were first introduced by Bate et al. (1995) in 
a SPH code. Sinks are massive particles that capture gas 
particles in their surroundings. They mimic the formation 
of unresolved compact objects, e.g. proto-stellar cores in the 
ISM, black holes in the ISM, central SMBHs in galaxies etc. 
Due to the very Lagrangian nature of the sink particle tech- 
nique, it has been extensively and exclusively used in SPH 
codes until Krumholz et al. (2004) extended its use to grid 
codes. The version in RAMSES (Teyssier 2002) is strongly in- 
spired by the Krumholz et al. (2004) numerical implementa- 
tion, and has already been presented in Dubois et al. (2010) 
and Teyssier et al. (2011), but we reproduce here the details 
of the numerical implementation to facilitate the discussion 
of our results. 



2.1.1 Seeding galaxies with BHs 

There are at least two scenarios for the formation of seed 
BHs. The first one invokes population III stars with zero 
metallicity. These can produce BH remnants as massive as 
10^-10^ Mq (Madau & Rees 2001; Heger & Woosley 2002; 
Schneider et al. 2002) that will eventually rapidly merge in 
their primordial halo to reach even larger masses. Another 
channel of BH formation is the direct collapse of matter 
from halos with very low angular momentum generating BHs 
as massive as lO'^ Mq (Loeb & Rasio 1994; Bromm & Loeb 
2003; Begelman et al. 2006). With the kpc-scales typically 



used in cosmological simulations of galaxy formation, it is 
pointless to try to follow the formation of these first seeds 
since this occurs on much smaller scales, but we can take 
these scenarios as guidelines for a sub-grid generation of seed 
BHs. 

BHs represented by sink particles are created in regions 
where the Jeans criterion is violated, i.e. in regions where the 
maximum level of refinement is reached and where the gas 
density is large enough to potentially produce a numerical 
instability, in other words where: 



Ax 



> Aj = 




(1) 



Here Aa; is the size of the smallest cell, Aj the Jeans length, 
Cs the sound speed and p the gas density. According to 
Truelove et al. (1997), the numerical stability of a gravita- 
tionally bound object is ensured if it is resolved with at least 
4 cells. With a mixed composition of matter (dark matter, 
gas, stars), Jeans stability is not trivial anymore, but we 
can reasonably assume that gas is the dominant source of 
gravitational potential inside dense collapsed objects, like 
galaxies, in our case. 

For numerical stability, each time that the Jeans cri- 
terion is violated we should spawn a sink particle with a 
mass corresponding to the depleted mass. However, in cos- 
mological simulations this leads to excessively large sink 
masses. The reason is that the gas is concentrated in galac- 
tic structures that are poorly resolved with kpc-scale res- 
olution. As a result an entire galactic disk can be defined 
by only a few Jeans-violating cells leading to excessively 
massive sink particles. To form sufficiently small seed BHs 
in the centres of the galaxies, we prefer to choose their 
initial mass, Mgcod, thereby introducing a free parameter. 
We set Msocd ~ 10^ Mq as the default value of our model 
in agreement with previous cosmological simulations (e.g. 
Booth & Schaye 2009). Despite choosing the seed mass, BHs 
are still spawned only in cells belonging to the maximum 
level of refinement and that verify equation 1. In section 4, 
the importance of the choice of the initial seed mass will 
be tested. One consequence of this self-controlled formation 
of the seed BHs is that they are not allowed to accrete gas 
when the Jeans criterion is violated. They can only accrete 
gas by a reasonable physical process such as Bondi accre- 
tion. With this prescription for initializing the mass of the 
seed BHs, it is conceivable that gas could be numerically 
violently Jeans unstable, but this issue is partially solved 
by the consumption of gas in the star forming process that 
temporarily restores gravitational stability. 

To avoid formation of sink particles in low density re- 
gions that are Jeans-unstable, we set a minimum threshold 
for the density p > po of gas that can create a new sink, 
where po is the same density threshold that we use for star 
formation. To make sure that sink BHs do not form before 
the very first stars form, we check that the local star den- 
sity p* calculated with a Cloud-in-Cell (CIC) interpolation 
verifies 



> 0.25, 



(2) 



before a new sink particle is spawned, where p is the local 
gas density. Note that these criteria are very similar to those 
employed by Bellovary et al. (2010), as they also confine the 
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formation of seed BHs to cold, dense, metal poor gaseous 
regions at the centre of galaxies. 

To obtain one BH per massive galaxy only, a halo finder 
is usually run on-the-fly during the simulation to check if 
candidate galaxies already host a BH (Di Matteo et al. 2005; 
Booth & Schaye 2009). We prefer a simpler, more direct, 
and computationally cheaper approach. To avoid creating 
multiple BHs inside the same galaxy, we ensure that each 
time a cell could potentially produce a sink particle (i.e. it 
verifies eq. (1)), it is farther than a minimum radius Vynin 
from all other pre-existing sink particles. This distance has 
to be larger than the typical size of galactic discs and smaller 
than the typical average inter-galactic distance. Test runs 
suggest that the choice rmin = 50 kpc produces very satis- 
factory results. 

In summary, a sink particle forms out of gas satisfying 
criteria on: Jeans instability, gas density threshold, stellar 
fraction threshold, and minimum distance from other BHs. 
Once the sink particle is created, it is split into several cloud 
particles with equal mass. Cloud particles are spread over a 
4:Ax radius sphere and positioned every 0.5 Ax in (x,y,z). 
The exact number of cloud particles in this configuration 
is therefore ndoud = 2109 per sink. These cloud particles 
are essentially created to probe the evolution of the region 
around the BH and provide spatially averaged quantities 
for the Bondi accretion formula. They move around on the 
finest time step scale (corresponding to the highest spatial 
resolution) and are destroyed and re-created around their 
sink particles with a given distribution at the beginning of 
every coarse time step (corresponding to lowest spatial res- 
olution). 

On the other hand, sink particles are only updated ev- 
ery coarse time step with quantities that have been evolved 
through the intermediate calculation with cloud particles. 
They are merged together (mimicking the BH merger) if 
they stand at a distance closer than AAx from each other. 
Mass is conserved in this process and momentum vectors of 
the old sink particles are simply added to compute the mo- 
mentum of the new sink particle. They are also the source 
AGN feedback. 

Finally we insist on the fact that BH positions and ve- 
locities are updated in the classical way used to update stan- 
dard particles such as DM particles, i.e. using the Particle- 
Mesh solver of RAMSES with CIC interpolation of particle 
masses into cells. No correction on their positions and ve- 
locities is done to force them to stay near their host galaxy 
(as could be done with a halo finder approach). Thus, weakly 
bound objects, such as BHs in galaxy satellites of large 
groups and clusters, may be stripped from their host galaxy. 
These BHs behave like star particles that tidal forces compel 
to populate the stellar halo of massive galaxies. 



(Bondi 1952) 



2.1.2 Accretion rate 

Since we do not resolve the accretion disks around BHs, 
whose sizes are sub-parsec even for the most massive ones 
(~ 10~^ pc according to Morgan et al. 2010 from micro- 
lensing estimates), we use the most common prescription 
that these BHs accrete gas at a Bondi-Hoyle-Lyttleton rate 



Mbh = 



47raG^M|HP 

(ci-h«2)3/2 



(3) 



where a is a dimensionless boost factor (q ^ 1), A/bh is 
the BH mass, p is the average gas density, Cs is the average 
sound speed, and u is the average gas velocity relative to 
the BH velocity. One of the major difficulties encountered 
with the computation of the relative gas velocity is that 
in cosmological runs, the ISM is poorly resolved and leads 
to galaxy scaleheights comparable to the resolution which 
is much larger than the scaleheights of galaxies in nature. 
Moreover, due to limited sampling of the gravitational force 
in the galactic disc, BHs can wander in their host galaxy. 
For this reason a BH close to the centre of a galaxy can feel 
the infalling material coming from the halo or the ISM at a 
relative velocity much higher than the typical velocity inside 
the bulge. Therefore because u is not a reliably measured 
quantity, we enforce the relative velocity to be no larger 
than an average gas velocity dispersion in the ISM which 
is assumed constant and equal to ttmax ~ lOkm.s"^ for our 
fiducial model (Dib et al. 2006). We will test the effect of 
varying this maximum allowed value on the properties of 
BHs. 

The average density p and sound speed Cs are computed 
around the BH using the cloud particles for this operation, as 
mentioned in the previous section. To compute the averages, 
the cell in which each cloud particle sits is assigned a weight 
given by a kernel function w, similar to the one used in 
Krumholz et al. (2004): 



w oc exp I 



(4) 



where r is the distance from the cloud particle to the sink 
particle and rK is the radius defined as 

( Ax/4 ren < Ax/A , 

rif = < ren Ax/4 < ren 5? 2Ax , 
[ 2 Ax r-BH > 2 Ax . 

The Bondi-Hoyle radius reH is given by: 

GA/bh 



(5) 



(6) 



where Cs is the exact sound speed in the cell where the sink 
lies. 

The accretion rate onto the sink is finally limited by its 
Eddington rate 



A^Edd = 



AnGMBHrrip 



(7) 



where ctt is the Thompson cross-section, c is the speed of 
light, rfip is the proton mass, and is the radiative efficiency, 
assumed to be equal to 0.1 for the Shakura & Sunyaev 
(1973) accretion onto a Schwarzschild BH. 

The accretion rate is computed at each time step and a 
fraction AfBHAf/ricioud of gas mass is depleted from the cell 
where the cloud particle lies and is added to that cloud par- 
ticle and to the sink particle. For each coarse time step of the 
simulation, cloud particles are re-scattered with equal-mass 
AfflH/ficioud. As the timestep does not depend on the accre- 
tion speed onto BHs and as low-density cells can be close 
to high density cells, a BH might remove more mass than is 
acceptable. To avoid dealing with negative or extremely low 
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gas densities and numerical instabilities arising from this, we 
do not allow any cloud particle to deplete more than 25% of 
the gas content in a cell. 

With large-scale cosmological simulations and the lim- 
ited typical kpc-scale resolution, we cannot resolve the scale 
and the dumpiness of the ISM that require parsec-like res- 
olution (Powell et al. 2011). To prevent the collapse of the 
gas from numerical instabilities and to take into account 
the mixing of the different phases in the ISM (cold and 
warm components), we use the polytropic EoS described 
in section 2.2. Applying this EoS means that it is impos- 
sible to track the 'true' density and sound speed in the 
ISM, thus the accretion rate onto BHs must be modified. 
Early works modelling the accretion rate onto BHs with 
such a polytropic EoS set the boost factor to a constant 
value a = 100 (Springel et al. 2005; Sijacki et al. 2007; 
Di Matteo et al. 2008). Here we follow the prescription from 
Booth & Schaye (2009) who show that a = [p/poY for 
p > po where po — O.lH.cm"^ is the threshold for star 
formation, and a — 1 for p ^ po give reasonable results 
compared to observational predictions. 

We insist on the fact that this polytropic EoS (equa- 
tion 16) has important consequences on the accretion rate 
onto BHs in high gas density regions: equation (3) turns into 
Mbh oc M|hP^^^, and the temperature dependence van- 
ishes. On the other hand, as soon as the cold gas component 
has been evaporated by star formation or feedback mecha- 
nisms giving p ^ Po in massive galaxies, the accretion rate 
of the BH is, by definition, the proper (a = 1) Bondi accre- 
tion rate. This a boost of the accretion rate is an artificial 
way of modeling the very fast accretion of gas within cold 
and gas-rich galaxies at early epochs, where the dumpiness 
of the ISM due to gas-disc fragmentation is unresolved in 
large-scale cosmological simulations. 



2.1.3 AGN feedback: Quasar and Radio mode 

It is believed that the feedback from AGN can proceed in 
two distinct modes. The quasar mode is essentially seen 
in the high redshift Universe and proceeds by emitting 
large amounts of radiation that can photo-ionize and heat 
gas. It is assumed in reionisation models of the IGM that 
quasars are an important contribution to the UV back- 
ground (Haardt & Madau 1996). The radio mode of AGN 
feedback, on the other hand, proceeds at lower redshifts in 
the cores of massive galaxy halos. The typical signatures of 
this radio mode are inflated cavities with strong magnetic 
fields and high levels of cosmic ray energy. 

Our aim is to treat self-consistently both modes in the 
simulation according to very simple prescriptions. It is be- 
lieved that the radio mode is preferentially triggered during 
low accretion-rate episodes, and that the quasar mode oc- 
curs when gas accretion takes place at rates comparable to 
the Eddington limit (Churazov et al. 2005; Merloni & Heinz 
2008). We use the ratio of accretion rate to its Eddington 
limit 

X=^ (8) 

MEdd 

as the criterion to determine which of the two AGN modes is 
active. Following Merloni & Heinz (2008), we take Xradio = 
10~^ as the value dividing the radio from the quasar mode. 



Above X > Xradio, the AGN undergoes quasar-like activ- 
ity with energy mostly emitted by photons. We model this 
mode by thermal injection of energy. Below x ^ Xradio, 
BHs smoothly accrete gas and provide a radio-mode feed- 
back which is modeled by our kinetic jet implementation. 
We point out that a similar approach has been taken 
by Sijacki et al. (2007), but they treat both modes as ther- 
mal inputs of energy with different injection radii. 

For both modes, we assume that a fraction ef of the 
radiated energy, Lr, is released to the ambient gas 

2 



-BaGN = ffir = eterMBHC 



(9) 



where ef is a free parameter that depends on the mode 
that is triggered by the accretion. As the energy is contin- 
uously released, time scales for dissipating energy by cool- 
ing can be sometimes far smaller than the hydro time step. 
This problem is often encountered in SN feedback mod- 
elling (Navarro & White 1993) and leads to a null dynam- 
ical impact on the surrounding gas. So as to impact the 
ambient medium, some authors release the AGN heating en- 
ergy only when a sufficient amount of gas has been accreted 
(Sijacki et al. 2007; Booth & Schaye 2009). Our modehng 
of the quasar mode [x > Xradio) as a heating mode is very 
similar to the approach adopted by Booth & Schaye (2009) 
(see also Teyssier et al. 2011): we store the rest mass en- 
ergy of gas accreted onto the BH until it would be enough 
to raise the temperature of the gas around the BH above 
T = lO^K. At that point we release the energy as thermal 
energy within a bubble of radius tagn around the BH with 
efficiency ef = 0.15. 

For the radio mode (x ^ Xradio) we model the AGN 
feedback with a jet-like outfiow with the same profile as 
in Omma et al. (2004) (see also Dubois et al. 2010). Meiss, 
momentum and energy are spread over a small cylinder of 
radius rj and height 2hj multiplied by a kernel window func- 
tion 



1 



exp 



r-2 



(10) 



where rcyi is the cylindrical radius, and where we impose 
r,j = hj = TAGN. The size of the jet, tagn, for the radio 
mode and the size of the bubble, tagn, for the quasar mode 
are parameter choices which we will test in section 4.1.6. 
Note that once the radius is chosen, it remains fixed for the 
duration of the simulation. By contrast, Sijacki et al. (2007) 
use prescriptions where the size of the radio bubbles depends 
on the amount of energy released and the gas density (see 
also Barai 2008). This different choice stems from the fact 
that these authors model the formation of bubbles by as- 
suming they are the result of large radio cocoons inflated by 
jets, whereas we attempt to directly model the jet. 
The mass deposition for the radio mode follows 

^ (r-cyi) „ 



Mj (r, 



cyl) 



\m 



-rjMBH ■ 



(11) 



where \\xp\\ is the integrated value of t/) over the whole cylin- 
der, and ri = 100 is an arbitrary value that represents the 
mass loading factor of the jet on unresolved scales. The value 
of rj adopted here corresponds to a sub-relativistic bipolar 
outflow with velocity 10 000 km/s, rather than that which 
is expected from a radio loud relativistic jet launched from 
the BH horizon. Thus our modelling should be interpreted 
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as a tentative description of the wind arising from a larger 
region surrounding the BH, and is expected to carry most of 
the momentum (see Omma et al 2004 for a more thorough 
discussion of this issue). Sucli a choice also allows one to 
keep the Courant time step of the simulation under control 
whist retaining a physically motivated model of the the jet 
outflow propagation on kpc scales. 

Mass is transferred from the central cell (where the BIf 
lies) to all the cells enclosed within the jet. Momentum, q, is 
deposited in outflowing opposite directions from the centre 
along the jet axis, according to 



|q.i|| (^-cyi) = Mj (rcyi) ||uj|| = ii!^MBHv/2^c. 



j.dr 



l|dr|| ' 
(12) 

where ||uj|| = (2efer/'7)"'^''^c is the velocity of the jet, (||uj|| ~ 
9487km. s^^ for ej — 1 and -q — 100), j is the spin vector of 
the BH and defines the jet axis, and dr is the distance vector 
from the centre of the BH. j is computed by adding the 
different contributions from the neighboring cells, sampled 
with the cloud particles, to the total angular momentum 



J = rriidri x ui , 



(13) 



where rrii and ui are the mass and velocity of the gas in 
the cell harbouring the cloud particle, i so that j — J/||J||. 
Finally the kinetic energy released within a single cell is 



Ej (rcyi) = 



2Mj (rcyi) ll^ll 



Eagt' 



(14) 



Integrating the energy over all the cells within the jet, we 
verify that the energy is strictly equal to i^AGN given in 
eq. 9. Energy efficiency et is a free parameter which we take 
to have different values depending on the AGN feedback 
mode, with et^ — 1 and ef,q = 0.15 our fiducial values for 
the radio and quasar mode respectively. 

High values of ef^r ~ 1 for the radio mode of AGN 
feedback are consistent with relativistic MHD simulations of 
BH accretion discs (De Villiers et al. 2005; Hawley & Krolik 
2006) for maximally spinning BHs. Using analytic argu- 
ments, Benson & Babul (2009) argue that balancing the 
spin up of BHs caused by gas accretion by angular momen- 
tum extraction through a jet, naturally leads to an equi- 
librium high spin value of the BHs, i.e. a/M = 0.93 corre- 
sponding to typical efficiencies tf.r ~ 0.1 and which seem to 
agree quite well with local observations (Allen et al. 2006). 
Of course such a picture is simplistic because of potentially 
rapid change of BHs spin during mergers, and the accretion 
mode of BHs (thin disks or ADAFs) is uncertain. However, 
it has the merit of yielding a very straightforward prediction 
of what should be the spin (and, thus the efficiency) of an 
isolated BH with an ADAF-like accretion disc, model which 
is consistent with our wind feedback assumptions. 

Unlike in Dubois et al. (2010), Dubois et al. (2011), in 
this implementation of kinetic AGN feedback (the radio 
mode), we allow for a time delay between the energy re- 
lease and its mass accretion similar to what Sijacki et al. 
(2007) do for thermal bubbles in their radio mode. The idea 
is that the energy is released into a kinetic jet when the 
BH has grown by more than AMd% of its mass. This pa- 
rameter, AMd, gives a relative, but artificial, control on the 
timescales of AGN feedback and their duty cycle, by allowing 



long periods during which the AGN is off and short periods 
during which the energy is released. 

Our approach for this dual radio/quasar mode of AGN 
feedback is obviously and voluntarily more simplistic than 
reality. As said before, the quasar mode involves the release 
of soft X-ray photons for which radiative transfer effects 
are not negligible. Also even during the quasar mode it is 
possible to get a faint radio detection, though it is not com- 
pletely clear whether such a signal would be coming from 
the remnant of a previously active radio mode or whether 
it is an intrinsic signal of the quasar mode. Jets are filled 
with a non-thermal cosmic ray component and strong mag- 
netic fields that could have important consequences on the 
dynamics of the jet. Moreover, the transition from the radio 
mode to the quasar mode very likely takes place at a differ- 
ent accretion ratio x than the transition from the quasar to 
the radio mode, reflecting the changing nature of the accre- 
tion disc. 



2.2 Modelling star formation and stellar feedback 

Gas in our simulation radiates energy by atomic collisions 
in a H/He primordial composition gas (Sutherland & Dopita 
1993) down to To = 10* K, so that it can collapse into DM 
potential wells to form galaxies. We also account for the en- 
hancement of cooling by metals released in SN explosions 
from massive stars. The metals are passively advected with 
the gas and a Solar composition of heavy elements is as- 
sumed. The minimum temperature To reached is not mod- 
ified by the presence of metals but they allow for a more 
efficient cooling. Heating from a UV background is consid- 
ered following Haardt & Madau (1996) during and after the 
redshift of reionisation which we take to be Zreion ~ 10.5. 

Star formation occurs in high-density regions with gas 
density p > po — O.lH.cm"^ using a random Poisson pro- 
cess to spawn star cluster particles, according to a Schmidt- 
Kennicutt law 



'iff 



(15) 



where p, is the star formation rate density, e* the star for- 
mation efficiency, and tg the gas fee-fall time. In these simu- 
lations we set the efficiency of star formation to e« = 0.02 in 
good agreement with observational surface density relation- 
ships of galaxies (Kennicutt 1998), and local giant molecular 
clouds (Krumholz & Tan 2007). Each star cluster particle 
has a mass of m« — poAx"^, reaching 3.6 10''' h^^.M0 for our 
most resolved simulation with Ax = 0.38 h^^. kpc. For nu- 
merical stability, we check that no more than 90% of the 
gas in a cell is depleted during the star formation process 
for numerical stability. 

We account for the mass and energy release from type 
II supernovae (SNe) assuming a Salpeter Initial Mass Func- 
tion (IMF). Using this IMF, 10% of the stars more mas- 
sive than 10 Mq end their life as type II SN releasing 10^^ 
erg per 10 Mq. Direct thermal input of energy from SN 
has been identified as an inefficient way of returning en- 
ergy back into the ISM because thermal energy is effi- 
ciently radiated away by gas cooling in high density regions 
(Navarro & White 1993). Approaches to circumvent this in- 
clude temporarily switching off gas cooling to allow the blast 
wave to propagate or kinetic energy feedback. The method 
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from Dubois & Teyssier (2008) implements the SN energy 
input by releasing mass, momentum, and kinetic energy lo- 
cally into the surrounding gas according to a Sedov blast 
wave solution. The explosion takes place 10 Myr after the 
birth of a star cluster particle and a fraction of the gas in 
the cell where the star particle resides is carried into the 
neighboring cells with a mass loading factor ~ 1- We as- 
sume that type II SNe release all their mass into the gas (no 
stellar remnant) with a y = 0.1 stellar yield, which is the 
fraction of primordial gas transformed into heavy elements 
and released back into the ISM. Our prescription does not 
take into account the energy and mass release from stellar 
winds (AGB stars), nor from long-lived type la SNe. 

In order to take into account the thermal impact of the 
heating of the ISM by SNe, we modify the temperature at 
high density p > po with a polytropic EoS 

T-n[^J-\ (16) 

where p is the polytropic index of the gas. The adopted 
value of p = 4/3 is comparable to the value obtained 
in Springel & Hernquist (2003) through analytic consider- 
ations of the multiphase structure of the ISM with stellar 
heating. This value of p = 4/3 does not rigorously ensure 
that gas will not fragment because of numerical instabili- 
ties (Truelove et al. 1997), as with this polytropic index the 
Jeans length is proportional to the gas density 

^■'-^°-K(nii^)"'''^P^- 

This last formula shows that at very high gas densities the 
Jeans length can be smaller than our minimum resolution 
and would cause spurious fragmentation of the gas. Fortu- 
nately, the gas cannot infinitely condense because of force 
resolution sampling and the star formation process that re- 
moves gas in the ISM. We did not choose the steeper, but 
safer, polytropic EoS, p = 2, because it leads to very thick 
galactic discs in simulations with kpc resolution. 



3 NUMERICAL ASPECTS 

We assume a flat ACDM cosmology with total matter 
(baryons-|-DM) density Qm = 0.26, baryon density Qi, = 
0.045, dark energy density Q,a = 0.74, fluctuation am- 
plitude at 8/i~^.Mpc (78 = 0.80 and Hubble constant 
Ho = 70km.s"\Mpc"^ consistent with WMAP 5-year data 
(Dunkley et al. 2009). We use several simulations with dif- 
ferent box sizes Lbox, number of initial DM particles A'dm, 
and minimum cell sizes Ax in order to test the convergence 
of our AGN feedback model with resolution. For a given 
Lbox size, we generate our most resolved initial conditions 
(ICs) and degrade them to obtain lower resolution ICs, so 
that for the same box size, the ICs will produce the same 
structures but with different numbers of DM particles. 

These simulations are run with the AMR code RAM- 
SES (Teyssier 2002). The evolution of the gas is followed 
with a second-order unsplit Godunov scheme for the Eu- 
ler equations. The Riemann solver for the flux computation 
at the cell interface uses a first-order MinMod Total Vari- 
ation Diminishing scheme to reconstruct the interpolated 
variables from their cell-centered values. Collisonless par- 
ticles (dark matter, stellar and sink particles) are evolved 



using a particle-mesh (PM) solver with a Cloud-In-Cell in- 
terpolation. 

Simulations refine the initial mesh by as many as 
7 levels of refinement, reaching a physical cell size of 
Ax = 0.38h"\kpc for our most resolved ICs and smallest 
box size Lbox ~ 12.5h~^.Mpc (simulations 256L12noAGN, 
256L12JH). Note that the fmax = 14 level of refinement is 
only reached at Oexp = (l-l-z)"^ = 0.8 for our most re- 
solved simulations. Because we enforce a nearly constant 
physical resolution (rather than a constant comoving res- 
olution), the highest refinement level triggered for a given 
redshift increases as the expansion factor grows with time, 
i.e. ^max 2 at ttcxp ~ 0.2, i'max 1 at ctexp — 0.4, etc. A 
cell is refined following a quasi-Lagrangian criterion: if more 
than 8 dark matter particles lie in a cell, or if the baryon 
mass exceeds 8 times the initial dark matter mass resolution. 
Lower thresholds for triggering refinement can be adopted 
to resolve the smallest halos with 10—1000 DM particles, us- 
ing sufficient force resolution (O'Shea et al. 2005). However, 
these can lead to excessive amplification of noise discreetness 
effects (Romeo et al. 2008). 

Our whole set of simulations using different box sizes, 
resolutions and sub-grid physics is summarized in table 1. 

Fig. 1 shows the gas density, temperature and metal- 
licity through a three-color composite image at two differ- 
ent redshifts for our most resolved simulation of a 50h~ 
Mpc box (256L50JH). with cooling, star formation, super- 
nova feedback and our fiducial model of AGN feedback, i.e. 
the dual quasar and radio modes. One can see the formation 
of hot and metal-rich bubbles at the intersection of filaments 
where galaxies collapse, already at high redshift z — 3, sug- 
gesting a pre-heating of halos before a cluster forms. The 
bubbles extend further into the IGM as the simulation pro- 
ceeds. 



4 CONSTRAINING THE AGN FEEDBACK 
MODEL 

In order to test the model for BH growth and its associated 
AGN feedback, we measure different statistically meaningful 
quantities and compare them to their observational equiv- 
alents. Quantities such as the cosmological density of BHs 
and the Magorrian et al. (1998) relationships that link BH 
masses to their host galaxy properties provide good con- 
straints on the coevolution of BHs and galaxies. We also 
test the effect of varying the AGN feedback parameters on 
the cosmic star formation rate (SFR per unit comoving vol- 
ume). 

The procedure adopted here is in many ways very sim- 
ilar to the one used in Booth & Schaye (2009) to constrain 
their AGN sub-grid model based on thermal energy inputs. 
As we share some common parameters with them, like the 
accretion rate boost factor a, and the thermal input of en- 
ergy for the quasar mode of AGN feedback, most of the 
parameter tests for the quasar mode will not be repeated, 
except for the AGN bubble size vagn because SPH and 
AMR codes treat their gas elements with a different numer- 
ical approach. Thus we set up our thermal (quasar) mode 
for AGN feedback to the values of Booth & Schaye (2009) 's 
best-fitting model: we assume an ef = 0.15 efficiency, which 
is also the value employed in Teyssier et al. (2011). In this 
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Table 1. Simulations performed with different sub-grid galactic models, different parameters for the AGN feedback mode, and different 
resolutions, (a) Name of the simulation, (b) Number of DM particles, (c) Mass resolution of a DM particle, (d) Size of the simulation 
box. (e) Minimum resolution reached at ^ = 0. (f) Presence of feedback from SNe. (g) Presence of AGN feedback: "BH" stands for the 
formation and growth of BHs without AGN feedback, "Jet" stands for the radio mode only, "Heat" stands for the quasar mode only, 
and "JET/HEAT" stands for the quasar and radio mode both triggered in the same simulation (see text (section 2.1.3) for details), (h) 
AGN feedback efficiency, (i) AGN energy delay, (j) Maximum relative velocity of the gas to the BH. (k) Mass loading factor of the jet. 
(1) Initial BH mass, (m) Size of the region for AGN energy input. 



(a) 


(b) 


(c) 


(d) 


(c) 


(f) 


(g) 


(h) 


(i) 


(j) 


(k) 


(1) 


(m) 


Name 




A/dm 


-^box 


Ax 


SN 


AGN 




AA/d 




V 


A/soed 


''AGN 






(Mo/h) 


(Mpc/h) 


(kpc/h) 








% 


(km/s) 




(M©) 




256L12noAGN 


256^ 


6.910'' 


12.5 


0.38 


Yes 


No 


- 


- 


- 


- 


- 


- 


256L12JH 


256^ 


6.9 10" 


12.5 


0.38 


Yes 


Jet/Heat 


1/0.15 


0/- 


10 


100/- 


10^ 


Ax 


64L25JH 


643 


3.5 10^ 


25 


3.04 


Yes 


Jet /Heat 


1/0.15 


0/- 


10 


100/- 


10^ 


Ax 


128L25noAGN 


128^ 


4.410* 


25 


1.52 


Yes 


No 


_ 


_ 


_ 


_ 


_ 


_ 


128L25BH 


128^ 


4.410* 


25 


1.52 


Yes 


BH 


- 


- 


10 


- 


10^ 


- 


128L25J 


128^ 


4.410® 


25 


1.52 


Yes 


Jet 


1 





10 


100 


10^ 


Ax 


128L25Je0.15 


128^ 


4.410® 


25 


1.52 


Yes 


Jet 


0.15 





10 


100 


10^ 


Ax 


128L25Je0.01 


128^ 


4.4 10* 


25 


1.52 


Yes 


Jet 


0.01 





10 


100 


10^ 


Ax 


128L25Jml 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 


1 


10 


100 


10^ 


Ax 


128L25JmlO 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 


10 


10 


100 


10^ 


Ax 


1 ooT OCT 1 nr\ 

lz8L25JvlUU 


128^ 


4.4 10* 


25 


1.52 


Yes 


Jet 


1 





100 


100 


lO'' 


Ax 


128L25Jvl000 


128^ 


4.4 10* 


25 


1.52 


Yes 


Jet 


1 





1000 


100 


10^ 


Ax 


128L25J»)10 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 





10 


10 


10^ 


Ax 


128L25J»;1000 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 





10 


1000 


10^ 


Ax 


128L25Js0.1 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 





10 


100 


10* 


Ax 


128L25JslO 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 





10 


100 


10^ 


Ax 


128L25J2dx 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 





10 


100 


10^ 


2Ax 


128L25J4dx 


128^ 


4.410* 


25 


1.52 


Yes 


Jet 


1 





10 


100 


10^ 


4Ax 


128L25H 


128^ 


4.410* 


25 


1.52 


Yes 


Heat 


0.15 




10 




10^ 


Ax 


128L25H2dx 


128^ 


4.410* 


25 


1.52 


Yes 


Heat 


0.15 




10 




10^ 


2Ax 


128L25H4dx 


128^ 


4.410* 


25 


1.52 


Yes 


Heat 


0.15 




10 




10^ 


4Ax 


128L25JH 


128^ 


4.410* 


25 


1.52 


Yes 


Jet/Heat 


1/0.15 


0/- 


10 


100/- 


10^ 


Ax 


256L25noSNAGN 


256^ 


5.5 lO'^ 


25 


0.76 


No 


No 














256L25noAGN 


256^ 


5.5 lO'^ 


25 


0.76 


Yes 


No 














256L25JH 


256^ 


5.5 lO'^ 


25 


0.76 


Yes 


Jet /Heat 


1/0.15 


0/- 


10 


100/- 


10^ 


Ax 


128L50noAGN 


128^ 


3.5 10^ 


50 


3.04 


Yes 


No 














128L50JH 


128^ 


3.5 10^ 


50 


3.04 


Yes 


Jet/Heat 


1/0.15 


0/- 


10 


100/- 


10^ 


Ax 


256L50noAGN 


256^ 


4.410* 


50 


1.52 


Yes 


No 














256L50JH 


256^ 


4.410* 


50 


1.52 


Yes 


Jet /Heat 


1/0.15 


0/- 


10 


100/- 


10^ 


Ax 



paper rather than focusing on the quasar mode, we concen- 
trate on testing our radio mode based on bipolar kinetic 
outflows, which is a modified version of the bipolar kinetic 
outflows in Dubois et al. (2011) and has never been modeled 
before in cosmological simulations of galaxy formation. 

To test the AGN sub-grid model, we fix the 'standard' 
galactic sub-grid models for star formation and SN feed- 
back (see section 2.2, e.g. star formation threshold po — 
O.lH.cm"'^, star formation efficiency e, = 0.02, Salpeter 
IMF, mass loading factor /u, = l, stellar yield y=0.1, poly- 
tropic index p=4/3) and we vary the parameters of the 
AGN feedback modeling (e/, 5Md, Umax, ?7, Msocd, rAGN, 
mode of AGN feedback). These tests are performed on a 
ibox = 25h~^.Mpc simulation box with 128^ DM particles 
(i.e. simulations with names starting with the prefix '128L25' 
in Table 1) that are sufficient to resolve halos with masses 
as small as ~ lO^^h~^.M0 (~100 particles) and as large 
as several ^ 10^"^ h~^.MQ. This choice of box size and res- 
olution is a good compromise between affordable computa- 



tional resources and resolution requirements. Finally conver- 
gence tests are performed with larger box sizes (128L50JH, 
256L50JH) and more (256L12JH, 256L25JH) and less re- 
solved simulations (64L25JH) for our radio/quasar mode 
with parameters from our best fitting model (simulation 
128L25JH, see table 1). 

In order to compare the simulated BH masses to their 
host galaxy bulge stellar mass as has been done for observa- 
tions, we must decompose the stellar surface density profiles 
into an inner bulge and an outer disc component. Our bulge- 
disc decomposition is detailed in Appendix A. 

Fig. 2 compares the BH comoving density as a func- 
tion of redshift while varying the parameters of the AGN 
feedback model. The value of the BH density in our local 
Universe from Shankar et al. (2004) is overplotted with its 
3(7 uncertainty. Any AGN feedback model that pretends to 
model the cosmological growth of BHs should remain close 
to this data point id z — 0. However, the latter does not en- 
sure that their growth is correct, as the observational point 
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Figure 1. Three color image of the 256L50JH simulation (see table 1) at 2 = 3 (upper panels) and z = (bottom panels) with a zoom 
on tlie largest halo (top and bottom right panels). Gas density is colour coded in magenta, gas temperature in cyan and gas metallicity 
in yellow. 



is only for redshift z — and the data at larger redshifts is 
more sparse. 

Fig. 3 shows the relationships between BH mass and 
their host galaxy stellar bulge mass as well as with the 
host's stellar velocity dispersion. These are good constraints 
for testing the coevolution of BHs and their galaxy mass 
content. We represent the average value of the distribution 
of stellar bulge masses for a given bin of BH mass. Obser- 
vational fits to BH mass (Mbh) versus stellar mass (Ms) 
(Haring & Rix 2004), and on BH mass (Mbh) versus stel- 
lar velocity dispersion (as) (Tremaine et al. 2002) are over- 
ploted with 3a uncertainties. Fig. 4 shows the cosmic SFR 



to test the impact of varying AGN feedback parameters on 
the history of mass assembly of galaxies. 



For these comparisons, we take simulation 128L25J as 
our reference model for the radio AGN feedback mode (see 
table 1), and each parameter of the AGN feedback model 
is varied one- by-one. Even though this procedure does not 
explore the full parameter space, and does not ensure that 
another set of parameters is possible, it allows us to test the 
validity of our fiducial model as well as the dependency of 
the results on the variation of its parameters. 



10 Y. Dubois et al. 





1 10 1 10 

l+z 1+z 




1+z 





(f) ■ 










L 128L25J 




128L25J2dx 




L 128L25J4dx 





10 



1+z 




1+z 1+z 



Figure 2. Comoving black hole mass density as a function of redshift. (a) Varying the efficiency ef for the jet mode, (b) Varying the 
AGN energy delay AM^j for the jet mode, (c) Varying the maximum relative velocity itmax i^or the jet mode, (d) Varying the mass loading 
factor for the jet mode, (e) Varying the BH initial mass Mgeed for the jet mode, (f) Varying the the AGN input size tagn for the jet 
mode, (g) Varying the the AGN input size t-agn for the heating mode, (h) Varying the mode of the AGN feedback. The dashed line is 
the average black hole mass density in our local Universe with its Sir uncertainty (grey shaded area) from Shankar et al. (2004). 
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og (l^sun) '09 ('<rn/s) log (Mj,^^) log ct^ (km/s) 






Figure 3. Panels (a) through (h) explore a variation of parameters in the same order as listed in the caption for fig. 2. In each panel 
we plot the black hole mass as a function of the stellar mass (left), or as a function of stellar velocity dispersion (right). Measurements 
are done at 2 = 0. We overplotted the observational laws as dashed lines from Haring & Rix (2004) for the Mbh-Mb relation and 
Tremaine et al. (2002) for the A/BH-fs relation with their 3cr uncertainties. The dotted line in the left-hand panel of (a) indicates the 
relationship between log Mbh and log Ms when Mbh = ^b- 
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l+z 1+z 

Figure 4. Panels (a) through (h) explore a variation of parameters in the same order as listed in the caption for fig. 2. Comoving SFR 
as a function of redshift. The dashed line corresponds to the 128L25noAGN simulation that does not include AGN feedback. 
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4.1 Parameter study 

4-1.1 Varying efficiency et 

We test the effect of varying the efhciency tf of the radio 
mode on the evolution of the cosmic BH density in fig. 2. (a). 
We compare three simulations with the implementation of 
the radio mode identical in every respect expect for the effi- 
ciency, ef which varies by approximately factors of ten from 
0.01 (128L25Je0.01) to 0.15 (128L25Je0.15) to 1 (128L25J). 
We also compare a run with black hole growth but no AGN 
feedback (128L25BH). Efficiency values lower than €{ = 1 
produce larger BH densities at z = than their observa- 
tional counterparts. The case for which no AGN feedback 
energy is released (128L25BH) produces what we take to be 
the maximum attainable BH density. Even small amounts 
of energy (ef = 0.15 and et — 0.01) prevent BHs from grow- 
ing to these maximum BH densities. The maximum possible 
efficiency ef = 1 predicts a slightly larger BH density than 
the data from Shankar et al. (2004), but is still within 3a 
error. 

The decrease in efficiency, e/ is compensated by 
larger accretion rates and more massive BHs, leading to 
the net result that the total amount of energy released 
by the AGN feedback is nearly independent of e/ (see 
also Booth & Schaye 2010). Fig. 5 substantiates this by 
showing the comoving cumulative AGN energy density as 
a function of redshift for different efficiencies. Indeed, even 
though AGN efficiency ef is allowed to vary by two orders 
of magnitude, the total amount of energy liberated at 2 = 
differs by less than a factor 2. This suggests that BHs ad- 
just their masses so that the total energy liberated can blow 
the gas out from the galaxies and stop the accretion of gas. 
We do not present the AGN energy density evolution for 
the other simulations presented in fig. 2 because aside from 
the simulations presented in panel (h), they all run with the 
same e/. Therefore the AGN energy density, eAGN, can be 
deduced from their pbh since eAGN = e/e^c^pBH, where we 
recall that the radiative efficiency er = 0.1 for all simulations 
and c is the speed of light. 

Decreasing the efficiency leads to more massive BHs, 
but, interestingly, the stellar masses of these galaxies and 
their stellar velocity dispersions are not significantly im- 
pacted by the efficiency (fig. 3. (a)). This is confirmed by the 
cosmic SFR seen in fig. 4. (a), which shows little difference in 
SFR for different values of the AGN feedback efficiency, es- 
pecially for ef = 1 and ef — 0.15. BHs regulate themselves as 
well as the gas content of their host galaxy by injecting the 
same quantity of energy, regardless of ef , to unbind the cold 
gas component. The small decrease in the SFR for ef = 0.01 
occurs because BH masses become comparable to their host 
galaxy masses so the BHs accrete gas instead of letting it 
form stars. This effect is more obvious when AGN feedback 
is not allowed but BH growth is permitted (128L25BH) : the 
SFR is suppressed by one order of magnitude because BHs 
are more massive than the entire stellar content of their host 
galaxy and they consume all the fresh gas available. 

We also plot in fig. 4, the cosmic SFR for the simula- 
tion (128L25noAGN) without AGN feedback nor BHs but 
including our standard sub-grid physics (cooling, star forma- 
tion, SN feedback). In this case, the SFR is systematically 
higher than any of the simulations including BH growth with 
or without AGN feedback, and the difference is clearer at 



60 




1 +z 

Figure 5. Cumulative comoving AGN energy density as a func- 
tion of redshift for different AGN feedback efHciencies ef. 

low and intermediate redshifts z = — 4. This shows that 
AGN feedback efficiently suppresses star formation in galax- 
ies, because it prevents gas overcooling and/or ejects large 
amounts of cold gas back into the Gircum- Galactic Medium 
(CGM). 

We remark that the jet velocity depends on the effi- 
ciency as itj (X s/TiJr], and that the cumulative momen- 
tum imparted by all BHs is Q oc ^/Ty/efeAON- Since eAGN 
is almost independent of ef (fig. 5), Q depends only on 
\fr\JTi. Thus, lower efficiencies produce higher total momen- 
tum providing a possible explanation for why self-regulation 
is weaker for lower efficiency, e/. However by varying the 
mass loading parameter, 77, in section 4.1.4, we show that 
self-regulation is controlled by the AGN feedback energy 
rather than momentum. 

This first set of experiments exploring variations in e/ 
suggests that AGN feedback is a necessary element for the 
self-regulation of the growth of BHs, and that a high value 
of the AGN feedback efficiency in the radio mode ef = 1 is 
required to obtain realistic results on the coevolution of BHs 
and galaxies. 

4.1.2 Varying AGN energy delay: AMd 

We allow for a time delay, represented by the AMd param- 
eter, before releasing AGN energy in the radio mode. This 
parameter prescribes that the BH must grow by more than 
a AMd fraction of its mass before releasing energy into a 
bipolar kinetic jet. We test three values for the time de- 
lay parameter for the radio mode: AMd ~ (128L25J), 
1% (128L25Jml), and 10% (128L25JmlO), At high redshift, 
varying the time delay has a negligible impact on the evo- 
lution of the BH density (see fig. 2.(b)) because BHs grow 
close to the Eddington accretion rate MEdd (see section 5). 
Thus, typical growth timescales of BHs are extremely short 
and even with a non-zero AMd, energy is released almost 
continuously. The BH growth timescale can be defined by 
tBH = Mbh/Mbh, and since MEdd oc Mbh, it simplifies to 

depending only on the Eddington accretion ratio x- a 
consequence, BHs accreting gas at high Eddington ratios 
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have a more continuous and hence immediate impact on the 
surrounding gas than BHs in a low accretion mode. 

However, at low redshift, results from simulations with 
different time delays diverge. When a delay is permitted 
(128L25Jml, 128L25JmlO), the final BH density &t z = 
is smaller than when energy is continuously deposited 
(128L25J). This is linked to easier BH self-regulation with 
larger AMd. Since BHs accrete gas at low redshift in a low- 
Eddington accretion regime, with non-zero AA/d there are 
larger stretches of time before a significant amount of en- 
ergy is released. As a result of this accumulation of energy, 
BHs release a comparable amount of energy to the AMd = 
case (128L25J) but in a shorter amount of time, allowing 
for fewer but stronger AGN luminosity bursts. Thus, BHs 
can more easily self-regulate with smaller duty cycles (Pope 
2011). 

The time delay parameterized by AMd also modifies the 
relationships between BH masses and their host galaxy prop- 
erties (fig. 3.(b)): BHs with masses Mbh > lO^Mo, sit in 
lower mass galaxies. This effect is stronger for AMd = 10% 
(128L25JmlO) for which stellar masses are reduced by an 
order of magnitude for the most massive galaxies compared 
to the simulation with a continuous (AMd = 0) injection 
rate (128L25J). The maximum stellar velocity dispersions 
of the host galaxies of the most massive BHs are reduced as 
a direct consequence of the reduction of the stellar mass in 
them. The effect is even more apparent for the intermediate 
BHs (10^ < Mbh < IO^Mq) in 128L25JmlO, where a clear 
deviation from the observational fit is observed. This sug- 
gests that the stellar velocity dispersion in massive galaxies 
is essentially controlled by the cold baryon content, and not 
by the total halo mass, which is hardly modified by the AGN 
feedback. 

Finally, in fig. 4.(b), we see that the intermediate and 
low redshift SFR depends a lot on the AMd parameter, 
which is not surprising since we saw that is also influences 
the Mbh-Ms relationships. Large AAfd more efficiently sup- 
presses the total SFR and has more impact on the gas con- 
tent in galaxies because they undergo shorter and stronger 
episodes of AGN feedback. This parameter study teaches us 
that not only is the amount of energy deposited important, 
but that the duration of the energy release plays a key role 
in unbinding the gas content of galaxies. Galactic gas ex- 
posed to a small deposit of energy can efficiently return to 
equilibrium by the gas dynamics and by the short cooling 
times involved in high density gas. 

4-1.3 Varying maximum relative velocity: Umax 

We measure the effect of varying the maximum allowed ve- 
locity itmax of the gas relative to the BH in the Bondi for- 
mula (equation 3). Increasing our fiducial value Umax = 10 
km/s (128L25J) up to Umax = 1000 km/s (128L25JvlQ00) 
decreases the overall BH densities (see fig. 2.(c)), but values 
are still consistent with the observations. As we would ex- 
pect, larger values of the relative velocity inhibit the growth 
of BHs. This effect already comes into play at high redshift, 
when mergers between galaxies are numerous sometimes re- 
sulting in violent excursions of BHs in their host galaxies, 
leading to large BH velocities relative to the dense gas com- 
ponent. 

This spurious effect comes from our inability to resolve 



the very small scales of the ISM within which BHs should 
be embedded. Some authors have circumvented this prob- 
lem by correcting the positions of BHs when they move too 
far from the gravitational potential well (private commu- 
nication with Volker Springel). Here we prefer to adopt a 
more straight-forward approach by limiting the maximum 
gas velocity relative to the BH in the formula for gas accre- 
tion (equation 3) rather than changing the position of the 
BH. This also allows us to follow the BHs that are ejected 
from their galaxies by strong tidal effects and gravitational 
friction, as material is stripped from galaxy satellites falling 
into massive halos. 

Fig. 3.(c) shows that different Umax produce dif- 
ferent BH masses and host galaxy stellar properties. 
When large relative velocities are permitted (128L25Jvl00, 
128L25vl000), galaxies tend to be more massive, and, as a 
result, have larger velocity dispersions. Even though the to- 
tal BH density at z = is hardly changed for different Umax, 
the SFR is very sensitive to Umax (fig. 4.(c)). Large Umax 
values tend to nullify the effect of the AGN feedback on the 
total SFR, and the SFR converges to the case where feed- 
back from AGN is not allowed. As a consequence Umax = 10 
km/s corresponds to a choice that allows for a non-spurious 
quenching of the accretion rate while keeping the dynamics 
of the BH particles completely self-consistent. 



4-1-4 Varying jet mass loading factor: rj 

The mass loading factor ry is a free parameter of the ra- 
dio mode for AGN feedback that controls the velocity the 
jet would have if it was propagating into a void. We com- 
pare three simulations with 77 = 10 (128L25J7)10), rj — 100 
(128L25J), and rj = 1000 (128L25J7?1000). As can be seen 
from fig. 2.(d), fig. 3.(d), and fig. 4.(d), the BH and galaxy 
properties are very insensitive to the adopted values of the 
mass loading factor of the jet. The reason is that the jet 
couples to the gas in its surroundings and AGN feedback 
becomes ineffective when the energy of the jet becomes com- 
parable to the binding energy of the gas. The only difference 
introduced by rj is that depending on its value, the jet will 
go more or less quickly into equilibrium with the gas, but 
as the liberated energy is the same regardless of rj, the same 
amount of gas is impacted by the jet. 

We insist on the fact that the total imparted momen- 
tum Q (X \/ 77/efeAGN loses its eagn dependence because the 
latter is constant given that BH densities are constant for 
different 77 (fig. 2.(d)) and ef is the same (ef=l) for the three 
simulations we are comparing. As a result, Q oc ^rj/ei as 
it was for the case with varying efficiencies (section 4.1.1). 
However since tf is a constant, momentum, Q, only depends 
on rj for the set of simulations compared in this section. 
Figs. 2.(d), 3.(d), and 4.(d) show that BHs and their host 
galaxy properties do not depend on rj, or equivalently on the 
momentum Q. Therefore we conclude that BHs and their 
host galaxy properties are only sensitive to jet energies, not 
their momenta. 



4-1.5 Varying initial BH mass: Afsocd 

We vary the BH initial seed mass by choosing values as 
smaU as lO"* Mq (128L25Js0.1) and as large as K^Mq 
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(128L25JslO). Simulation 128L25J has Meoed = 10^ Mq. As 
seen on fig. 2.(e), this has little effect on the final (z — 0) 
BH density. Differences appear at high redshift, when most 
of the contribution to the BH density comes from Bffs with 
mass close to their initial seed mass. Thus, the choice for the 
initial seed BH mass has an important impact on the BH 
density at high redshift but memory of it is rapidly erased 
when BHs grow to values larger than their initial mass. 

The Magorrian relationships are almost unchanged by 
different choices for Msocd (fig- 3.(e)). But a careful inspec- 
tion of the cosmic SFR (fig. 4.(e)) shows that the SFR is 
slightly different at high redshift 2 = 2 — 4: lower seed BH 
masses result in a larger SFR because it takes more time 
for BHs to reach a self-regulated equilibrium with their en- 
vironment. 

4-1.6 Varying AGN input size: tagn 

We vary the size of the AGN energy input region for both 
the radio and the quasar mode to deduce its effect on BH 
self-regulated growth. 

First we explore the impact of varying the extent of 
the jet in the radio mode through a series of three simula- 
tions where the AGN input region varies from tagn = Aa; 
(128L25J), to TAGN = 2 Ax (128L25J2dx), to tagn = 4Aa; 
(128L25J4dx). Fig. 2.(f) shows that large jet sizes for the 
radio mode lead to large BH densities. Since the energy is 
spread over larger regions when we increase the jet size, 
the gas close to the BH is impacted less by larger jets, 
thereby accreting more easily onto the BH, hindering its 
self-regulation, and producing larger BH densities. 

Fig. 3.(f) shows that jets with sizes larger than Ax lead 
to larger BH masses and lower host galaxy stellar masses, 
with unrealistic values compared to observations. Also, the 
SFR is strongly suppressed when larger jet sizes are chosen 
(fig. 4.(f)). Thus, large energy injection regions for the radio 
AGN feedback mode have more impact on galaxy formation 
because BHs become more massive and as a result inject 
more energy to the surrounding gas. However, as the BHs 
are much more massive than what is predicted from the 
Mbh-Ms relationship, they must release more energy to self- 
regulate. 

For the quasar mode the behavior is different. Again, 
we ran three simulations to explore the effect of chang- 
ing the size of the energy input region: from tagn = Aa; 
(128L25H), to TAGN = 2 Ax (128L25H2dx), to tagn = 4Aa; 
(128L25H4dx). Fig. 2.(g) shows that the trend of BH den- 
sity with bubble size is non-linear. Doubling the radius 
(rAGN = 2Aa;) gives a similar BH density evolution as ob- 
tained with tagn = Ax down to redshift z = 1.5 but shows 
a drop in BH density below this redshift. Increasing the ra- 
dius by four (rAGN = 4Aa:), on the other hand, leads to a 
larger BH density at high redshift, which converges to the 
same value as the tagn = Aa; case a.t z = Q. 

For the Magorrian relations, choosing tagn = 2Ax 
rather than tagn = Ax for the quasar mode leads to smaller 
stellar velocity dispersions but very similar stellar masses. 
It seems that for the tagn = 2 Ax case, even though the 
SFR is significantly decreased compared to the tagn = Ax 
case (see fig. 4.(g)), the decrease in the BH density keeps the 
BH mass versus host galaxy stellar mass relation unchanged. 
However a larger energy injection region (rAGN = 4Ax) has 



a dramatic impact on the final galaxy stellar masses and 
the evolution of the SFR. Both are significantly diminished. 
Similar behavior has been found by Booth & Schaye (2009), 
where increasing the number of neighboring SPH particles 
affected by the AGN bubble decreases the SFR and increases 
the BH density. 

Both the radio and quasar modes experience a decline in 
SFR as the size of the injection region (figs. 4.(f) and (g)) in- 
creases because large energy injection regions extend to less 
dense regions which are easier to impact. By blowing out the 
reservoir of hot gas, accretion onto galaxies and hence SFR is 
suppressed. However the self-regulation of BHs is somehow 
very different for the two modes. The difference resides in 
the very nature of energy deposit. Jets put momentum and 
kinetic energy into the gas and eventually some of this en- 
ergy is transformed into thermal energy through shocks, but 
the more extended the jet, the weaker the shock. Fig. 6 il- 
lustrates this effect in the high redshift Universe; the quasar 
mode with large tagn = 4 Ax (128L25H4dx) inflates larger 
and hotter bubbles than the radio mode with the same initial 
jet extent (128L25J4dx). As the accretion rate is very sensi- 
tive to the temperature of the gas Mbh oc T~^'^ , it is more 
difficult for jets than for thermal bubbles to self-regulate the 
growth of BHs. This is why BH densities for the jet mode 
with rAGN = 4 Ax are larger than for the heating mode. 

As a final remark, these particular numerical exper- 
iments demonstrate that the injection of energy through 
AGN feedback (but it is true for any type of feedback, 
see Dalla Vecchia & Schaye 2008 for a similar discussion 
about SN feedback) is a delicate process that cannot be 
naively decoupled from the gas dynamics up to large dis- 
tances, and must be handled with great care. 

4-1-7 Comparing radio mode and quasar mode 

We compare the choice of using the radio mode (jet mode) 
to the quasar mode (heating mode), as well as to a combi- 
nation of both modes. We found a set of parameters (e/ = 1 
(radio) 0.15 (quasar); AMd = 0; «max = 10 km/s; 77 = 100; 
Msocd ~ 10^ Mq; rAGN = Ax) consistent with observations 
for the radio mode and quasar mode used individually and 
use this same set of parameters for the dual radio/quasar 
mode (128L25JH). With the dual radio/quasar mode of 
AGN feedback, the feedback from BHs switches to the radio 
mode when the Eddington accretion ratio x ^ Xradio and to 
the quasar mode when x > Xiadio- 

Fig. 2.(h) shows that the radio (jet) mode (128L25J) 
produces a larger BH density than the quasar (heating) 
mode (128L2H), suggesting that the quasar (heating) mode 
is slightly more efficient at self-regulating BH growth (al- 
though a slightly smaller AGN efficiency, £/, of the quasar 
mode would reduce this difference). As expected, the dual 
mode feedback (128L25JH) gives a BH density value which 
is between the value for individual modes. 

There are some 'jumps' in the evolution of the BH den- 
sity at z = 4 and z = 1.5 for the radio (jet) mode (fig. 2.(h)), 
which are not present for the quasar (heating) mode. This 
effect is also seen in the SFR evolution (fig. 4.(h)), both for 
the radio (jet) and quasar (heating) mode. It comes from the 
triggering of a new maximum level of refinement at these 
redshifts. With one more level of refinement, the force on 
small scales is better resolved and the gravitational poten- 
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Figure 6. Projected temperatures of the 128L25JH (left), 128L25J4dx (middle), 128L25H4dx (right) simulations at 2 = 3. The color 
code is in log(K) units. Images are 6.25 h~^.Mpc physical on a side. 



tial well is deeper. The gas can therefore get compressed to 
higher densities, leading to significant boosts in its star for- 
mation rate and the accretion rate onto the BHs. This effect 
is less pronounced for the quasar (heating) mode where the 
heating supplies an additional pressure support to the gas, 
tending to erase this spurious numerical feature. 

Fig. 3.(h) shows that the stellar mass in the quasar 
(heating) mode (128L25H) is slightly reduced compared to 
the radio (jet) mode (128L25J), which is another signature 
of a slightly more efficient feedback in the quasar (heating) 
mode. The mode of AGN feedback has very little effect on 
the BH mass versus stellar velocity dispersion. The com- 
bination of the two modes show deviations, especially at 
low Ma and (Tb, from observed BH mass versus stellar mass 
relationships and observed BH mass versus stellar velocity 
dispersion relations. However in fig. 9 we show that this is a 
consequence of the limited numerical resolution, and that we 
converge to the observational measurements with increased 
resolution. 

Fig. 4.(h) shows that the SFR for the different AGN 
feedback modes are almost undistinguishable. A small dif- 
ference can be seen at low redshift 2 = — 1 where the quasar 
(heating) mode seems to more efficiently suppress the total 
SFR than the radio (jet) mode. The latter explains the dif- 
ference seen in the BH mass versus stellar mass relationships 
(fig. 3.(h)). 

Finally, this parameter study has allowed us to choose 
the best fitting parameters for our dual radio/quasar AGN 
feedback model compared to observations, namely the pa- 
rameters used in the 128L25JH simulation (see table 1). 



4.2 Resolution study 

In order to test the convergence of our models, we vary the 
resolution, by changing the DM mass, cell size, and box size. 
We run five different simulations with our fiducial model for 
AGN feedback (256L12JH, 256L25JH, 128L25JH, 64L25JH, 
256L50JH, 128L50JH), with three different box sizes Lbox = 
12.5h"^Mpc, Lbox = 25h~\Mpc and Lbox = 50h"\Mpc, 
and four different resolutions {Mdm ~ S.SIO^M©, Aa; = 
3.04h-\kpc} (64L25JH, 128L50JH), {Mdm = 4.41O*M0, 
Ax- = 1.52h"^kpc} (128L25JH, 256L50JH), {A/dm = 
5.5 lO'^ Mq, As: = 0.76h-^kpc} (256L25JH), and {Mdm = 
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Figure 7. Comoving black hole mass density as a function of 
redshift for different box sizes and resolutions. The grey shaded 
area is the black hole mass density in our local Universe with its 
3cr uncertainty from Shankar et al. (2004). 



6.910® Mq, Ax- = 0.38h"\kpc} (256L12JH) (see table 1 for 
details) . 

The BH densities shown in fig. 7 slowly converge to the 
same value at 2: = when the resolution is increased. Even 
though low resolution simulations (64L25JH and 128L50JH) 
are within the observational 3 o error bars, they tend to 
underestimate the BH density at all redshifts compared to 
more resolved simulations. Intermediate resolution simula- 
tions with {Mdm = 4.410® M©, Aa; = 1.52h"\kpc}, which 
correspond to runs 128L25JH and 256L50JH, have already 
converged at 2 = and differ only by ~ 10% at z=0 from 
the simulation (256L25JH) with one additional refinement 
level. However, at high redshift, the difference is larger be- 
cause galaxies in the field with intermediate BH masses con- 
tribute more to the total BH density (fig. 8), and some of 
these galaxies are not resolved in simulations 128L25JH and 
256L50JH with lower resolution. 

The same convergence can be seen in the BH mass 
versus stellar mass relationships (fig. 9, left panel). When 
the resolution is increased, these relationships converge to 
the same value close to the observations from Haring & Rix 
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Figure 8. Comoving black hole mass density as a function of 
redshift with contributions from different BH mass ranges. Upper 
panel is for the 256L25JH simulation, and bottom panel is for 
the 256L50JH simulation. The grey shaded area is the black hole 
mass density in our local Universe with its Scr uncertainty from 
Shankar et al. (2004). 
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Figure 9. Black hole mass as a function of the stellar mass (left 
panel), or as a function of stellar velocity dispersion (right panel) 
at z = for different box sizes and resolutions. The color code is 
the same as for fig. 7. We overplotted the observational laws as 
dashed lines from Haring & Rix (2004) (left) and Tremaine et al. 
(2002) (right) with their Scr uncertainties. 

(2004). There is a departure from the observational con- 
straint for the least massive galaxies that reside in halos 
that are barely resolved. For PM codes, DM halos with more 
than ~ 1000 DM particles are followed with sufficient force 
resolution (O'Shea et al. 2005; Heitmann et al. 2008). Thus, 
the break observed at the low galaxy mass end corresponds 
to this low limit. 
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Figure 10. Comoving SFR as a function of redshift for different 
box sizes and resolutions. The color code is the same as for fig. 7. 
The light grey circles with error bars correspond to observational 
points from Hopkins & Beacom (2006). 



For example, the 256L25JH simulation has a DM mass 
resolution A4-dm ~ 5.5 10^ h~^.M0, which gives a minimum 
DM halo mass Afh.min ~ 8 10^° M© that corresponds to a 
total gas content of /6Mh,miii ~ 10^° Mq. A non-negligible 
fraction of the baryon content is locked into stars. Assuming 
25 % for this fraction, we obtain A/s,min ^ 2.51O^M0, i.e. 
the value of the bulge stellar mass where the break appears 
in the Mbh-Ms relationship. Indeed, the same behavior is 
seen when relating the BH mass to the stellar velocity disper- 
sion of their host galaxy: low-mass BHs hosted by galaxies 
with low stellar velocity dispersion show a significant devi- 
ation from the constraints of Tremaine et al. (2002) (fig. 9, 
right panel). 

Fig. 10 shows the cosmic SFR for simulations with our 
fiducial model for AGN feedback and for different resolu- 
tions and box sizes. It appears that this quantity is strongly 
dependant on DM mass resolution at high redshift because 
of two effects. When increasing the resolution, smaller ha- 
los which contribute a lot to the total star formation rate at 
high redshift are resolved. Furthermore, structures that were 
already resolved at lower resolution collapse earlier when 
resolution is increased, thus, forming stars earlier. This ef- 
fect is well constrained and, if one has sufficient resolution, 
the numerical solution converges to its analytical prediction 
(Rasera & Teyssier 2006). 

We point out that Tree codes or Tree-PM codes, such 
as the GADGET code (Springel 2005), have a better force 
resolution in the early Universe than PM codes. The former 
therefore can follow the formation of smaller halos (down to 
10 DM particles) with equivalent initial conditions. This is 
why, in general, the convergence in the SFR with such codes 
is more rapidly obtained, even though 10 particles per halo 
are not sufficient to properly treat the gas dynamics. 

Even though the convergence for a given box size is not 
reached, the SFR for the box size 25h~^.Mpc slowly con- 
verges at z = 0. It increases by a factor 4 from 64L25JH to 
128L25JH, and by a factor 2 from 128L25JH to 256L25JH. 
Changing box size at constant resolution has also some non- 
negligible effect on the SFR at low redshift. The reeison is 
that at low redshift the SFR is essentially dominated by 
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Figure 11. Comoving SFR as a function of redshift for differ- 
ent box sizes Lbox = 12.5h~^.Mpc (red), Lbox = 25h~^.Mpc 
(black), Lbox = 50li-l.Mpc (blue) including SN and AGN 
feedback (solid), including SN feedback and no AGN feedback 
(dashed), and without SN or AGN feedback (dotted). The light 
grey circles with error bars corresponds to observational points 
from Hopkins & Beacom (2006). 
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massive galaxies. Thus, with larger box sizes, very massive 
clusters of galaxies are more likely to be present (Davis et al. 
2011). Finally, we conclude that getting convergence for the 
SFR is extremely difficult, because it requires both a large 
box size and a small DM mass resolution, which can only be 
achieved with tremendous computational power. 

It is interesting to note that the DM resolution effect on 
SFR is also present for the simulations without AGN feed- 
back (fig. 11), which proves that this effect is uncorrelated 
to AGN feedback but only to DM mass resolution and box 
size. AGN feedback is most efflcient at suppressing the SFR 
at low redshift when massive structures such as groups and 
clusters of galaxies are formed. At high redshift the SFR is 
dominated by galaxies in the field which are not progenitors 
of groups or clusters of galaxies, thus, the effect of AGN 
feedback is less visible. 

There is an interesting behavior of the simula- 
tion without AGN feedback and without SN feedback 
(256L25noSNAGN, dotted line in fig. 11) compared to 
the simulation without AGN feedback but with SN feed- 
back (256L25noAGN, black dashed line in fig. 11). At high 
redshift SN feedback reduces the SFR, because galaxies 
form large-scale galactic winds (Springel & Hernquist 2003; 
Dubois & Teyssier 2008) that remove some baryons from 
them and prevent some gas from collapsing into them. But as 
time goes by, structures become more massive and the ram- 
pressure confinement from the halo becomes higher, pre- 
venting galactic winds from escaping the discs. As a result 
they develop galactic fountains (Dubois & Teyssier 2008). 
The SN feedback also enriches the gas with metals and en- 
hances the gas cooling rates, developing stronger accretion 
flows and SFRs at late times, when metals are confined in 
halos (Dubois et al. 2011). Thus the feedback of SN has a 
negative impact on SFR at high redshift but a positive effect 
at low redshift due to metal enrichment. 



Figure 12. For each plot: black hole mass as a function of the 
stellar mass (left panel), or as a function of stellar velocity disper- 
sion (right panel). Measurements are done at different rodshifts 
as labelled in the upper right panel. Top plot corresponds to the 
256L25JH simulation, and bottom plot to the 256L50JH simula- 
tion. We overplotted the observational laws as dashed lines from 
Trcmaine et al. (2002) and Hiiring & Rix (2004) with their 3(7 
uncortainties. 



5 REDSHIFT EVOLUTION OF BH 
PROPERTIES 

The way BHs acquire their mass and how they liberate 
energy to the gas is important for constraining their co- 
evolution with their host galaxy. Observations of the rela- 
tionships between BH masses and galaocy masses at high 
redshift are extremely difficult because the luminosity of 
host galaxies of massive BHs is dominated by the AGN 
component. Despite these difficulties, an increasing amount 
of data is suggesting that the Mbh/A^s relationship shows 
some positive evolution with redshift (McLure et al. 2006; 
Peng et al. 2006; Shields et al. 2006; Salviander et al. 2007; 
Bennert et al. 2010; Decarh et al. 2010; Merloni 2010). Sim- 
ulations can provide insight on the time evolution of these 
relationships. 

Fig. 12 shows the BH mass versus bulge stellar mass at 
different redshifts for the 25 h"\Mpc (256L25JH) and the 50 
h~\Mpc (256L50JH) simulations. At high redshift, the BHs 
are more massive than the z = Q Magorrian relation would 
predict given the stellar mass of their host galaxy. This is 
supported by observations (see Merloni 2010). We evaluate 
this deviation from the z = relationships by measuring the 
median of the distribution of the AfaH/Mg measurements for 
the 256L25JH and the 256L50JH simulations in fig. 13. We 
explore the effect of removing from the sample BHs with 
masses smaller than a mass threshold to see if there is a BH 
mass for which the deviation is most pronounced. Like the 
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Figure 13. Median value of the BH mass (A/bh) over its host 
stellar bulge mass (Ms) as a function of redshift for the 256L25JH 
simulation (upper panel) and the 256L50JH simulation (bottom 
panel). Different colors correspond to different BH mass cut-offs 
used to get the median value of Mbh/Mb- The vertical error bars 
correspond to first and third quartile limits of the distribution 
of points. The dotted line is the trend from a fit to observa- 
tions in Decarli et al. (2010) and the dashed line is the trend 
from Merloni (2010). 



observations (Decarli et al. 2010; Merloni 2010) for which 
BH masses are larger than a few 10^ M© , we observe a posi- 
tive trend with redshift of the Mbh/Ms ratio. In our simula- 
tions, the ratio is larger for more massive BHs, but the trend 
is independent of the BH mass threshold. Quantifying the 
trend, we find Mbh/Ms <x {1 + z)"= with = 0.42 ± 0.09 
for the 256L25JH simulation and a, = 0.42 ± 0.06 for the 
256L50JH simulation when fitting our simulation data on 
BHs with masses larger than > 5 10^ M0 between redshifts 
a = to z = 3 . This is in relatively good agreement with the 
value Qs = 0.68 ± 0.12 measured in the observational data 
by Merloni (2010), and with numerical simulations from 
Di Matteo et al. (2008) [a^ = 0.5) and Booth & Schaye 
(2011) {as =0.52 ±0.05). 

The increase in Mbh/Ms reflects the different accretion 
modes onto the BH and the gas content and properties at 
different redshifts. Fig. 14 illustrates how massive BHs grow 
through time, with very fast accretion of gas at high redshift 
due to the presence of a cold and dense ISM. The accretion 
proceeds by bursts accompanied by large releases of AGN 
energy that temporarily delay the accretion onto the BH. At 
high redshift, two kinds of accretion occur: the accretion of 
a diffuse component that can eventually shock and virializc 
the gas in the halo, as well as accretion of dense filaments 



Figure 14. BH mass as a function of redshift for the most massive 
BH in the 256L12JH simulation (upper panel) and the logarithm 
of the ratio of its accretion rate to the Eddington accretion limit 
(bottom panel). 

of gas (Keres et al. 2005; Ocvirk et al. 2008; Brooks et al. 
2009; Dekel et al. 2009). The gaseous filaments feed galaxies 
so that their BHs can grow to larger masses, pre-heat their 
proto-cluster environment and remove gas thereby halting 
the adiabatic contraction in the proto-cluster cores. 

For simulations 256L50JH and 256L50noAGN, fig. 15 
shows the fraction of baryons in galaxies in the form of a cold 
gas component with gas density larger than > O.lH.cm^"^ 
for different stellar masses at different redshifts. It appears 
that more massive galaxies have lower gas fractions that de- 
cline with time. This can be explained by two effects. Firstly, 
galaxies efficiently consume their gas to form stars without 
replenishing their cold gas content quickly enough to main- 
tain a constant specific star formation rate. Secondly, AGN 
feedback reduces the amount of cold gas available in galax- 
ies by ejecting dense material into the CGM. Hence AGN 
feedback coupled to a vigorous consumption of gas via star 
formation reduces the gas content in galaxies, resulting in a 
much lower accretion rate at low-redshift, where BHs enter 
a low-Eddington accretion regime (fig. 14). 

The behavior of the accretion rate shown in fig. 14 for 
a single BH is common to BHs of a large range of masses. 
We represent the number-weighted diagram of the BH ac- 
cretion luminosity Lace = Mbhc^ versus the BH mass at 
different redshifts in fig. 16 for the 50 h~^.Mpc simulation 
(256L50JH). At high redshift {z = 4), the accretion pro- 
ceeds in a high accretion regime, where most BHs accrete 
very close to their Eddington accretion rate and release their 
energy in a 'quasar' mode. Some of the most massive BHs 
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Figure 15. Average gas mass in the disc (Mism) over the total 
mass of gas plus stars (Mism + Ms) as a function of the total stel- 
lar mass (Ms) for different redshifts z = 4, 2: = 3, ^ = 2, z = 1, 
z = 0.5, and z = from top to bottom for the 256L50JH simu- 
lation (solid lines) and for the 256L50noAGN simulation (dashed 
lines). 



have already entered a low-accretion regime suggesting that 
the most massive objects are the first to self-regulate their 
gas content. Then, as the simulation evolves, more and more 
BHs enter the low- accretion regime providing a 'radio' mode 
of AGN feedback. The core of the distribution which is lo- 
cated at ~ 10~^-1 MBdd at z = 4 is at ~ 10~^ Msdd at z = 2, 
~ 10~^ Mfidd at 2 = 1, and ~ 10~* A^Edd at z = 0, with very 
fewer and fewer Eddington- limited BHs at lower redshifts. 
We observe a lower-limit trend that goes like Mgjj produced 
by a combination of the minimum density and maximum 
temperature reached in the CGM very close to the galaxy 
hosting the BH. This lower bound evolves with time be- 
cause of both a rarefaction of the gas and an increase in 
temperature as haloes get more shock-heated as they grow 
in mass (Birnboim & Dekel 2003; Dekel & Birnboim 2006; 
Birnboim et al. 2007). 

These diagrams suggest that the feedback from AGN 
at high redshifts is essentially dominated by a quasar mode 
(x > 10"'^), whereas at low redshift a radio mode (x 10~^) 
prevails within the core of massive structures (Dubois et al. 
2010). 

Fig. 17 shows the AGN luminosities Lagn = 
e/erMBHC^ for different redshifts and two different box sizes 
25 h-\Mpc and 50 h-\Mpc (256L25JH, 256L50JH). Recall 
that for the dual quasar /radio AGN mode, e/ depends on 
the accretion rate to Eddington ratio x- The bright end of 
the high redshift (2 = 4) AGN luminosity function is domi- 
nated by the quasar mode because most of the BHs accrete 
gas at a high Eddington rate (see fig. 16). At intermediate 
redshifts, z = 2 and 2 = 1, the bright end of the AGN lumi- 
nosity is marginally dominated by the quasar mode, but the 
transition from the quasar mode to the radio mode appears 
at larger luminosities. Emitting the largest amounts of en- 
ergy, massive BHs start to strongly deplete the gas content 
of their host galaxies and enter a radio mode AGN regime. 
At 2 = most of the BHs have reached a very quiescent 
phase for gas accretion and quasar mode feedback is almost 
imperceptible. The last remaining AGN quasars are interme- 



diate mass BHs (see fig. 16), with most of the very massive 
BHs in a radio mode. 

It is worth noting that as a results of DM mass resolu- 
tion, the AGN luminosity functions show extrema. However 
the slope of the bright-end seems relatively independent of 
the resolution. 



6 DISCUSSION AND CONCLUSIONS 

We have designed a new self-consistent model of BH growth 
and AGN feedback for hydrodynamical cosmological simu- 
lations with a dual jet/heating mechanism that accounts for 
the radio mode and the quasar mode. BHs are seeded at an 
early stage of the formation of galaxies. They grow by suc- 
cessive mergers and by accretion of gas at a Bondi (1952) 
rate. Some of the rest-mass accreted energy is converted into 
energy for the AGN feedback mechanism. At high accretion 
rates defined by the Eddington ratio x > Xradio, BHs enter 
a quasar mode, with energy liberated as thermal energy. At 
lower accretion rates x ^ Xradio, a radio mode is triggered 
with injection of mass, momentum and kinetic energy within 
a bipolar jet. 

The parameters of the radio mode for the AGN feed- 
back model are tested to reproduce the observations of the 
2 = cosmic BH density and the Magorrian et al. (1998) 
relationships for A^bh-Ms and Mbr-cts- We find the follow- 
ing behavior of BH growth upon varying the parameters of 
the radio AGN feedback mode: 

- AGN feedback efficiencies lower than tf = 1 lead to 
larger and unrealistic BH masses that overshoot the ob- 
servational predictions for the cosmic BH density and the 
Magorrian relationships. BHs grow to larger masses to in- 
ject the same total amount of energy in order to self-regulate 
their growth. 

- Introducing a time delay in the AGN feedback, as op- 
posed to an instantaneous energy deposit, increases the ef- 
fectiveness of AGN feedback effect on the gas and the growth 
of BHs, as a more important energy release occurs over a 
shorter period of time. However, this increased efficiency 
tends to exaggerate the effect of AGN feedback, i.e. it un- 
derestimates the BH density and overestimates the Mbh/Ms 
ratios. 

- The mass loading factor of the jet, i.e. the velocity of 
the jet, has a negligible impact on the results. 

- The choice of the seed BH mass is relatively unimpor- 
tant except at high redshift when BH masses are compara- 
ble to the total stellar mass of galaxies. But, the choice of 
the seed mass is quickly forgotten as BHs self-regulate their 
growth. 

- The size of the jet for the radio mode, as well as the size 
of bubbles for the quasar mode, must be chosen carefully. 
Large jets and bubbles deposit energy far away from the 
galaxy and have a harder time impacting the dense gas, 
and, thus, self-regulating the growth of BHs. We find that 
the size of the energy injection region must be as close as 
possible to the minimum physical resolution of the code. 

- We find that the radio mode alone requires a larger 
energy efficiency ef,r = 1 than the quasar mode alone 
ef,q = 0.15 to match the data from observations. The ther- 
mal mode couples more efficiently to the gas than the kinetic 
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Figure 16. Number-weighted hiistogram as a function of the BH mass and their accretion luminosity for the 256L50JH simulation. The 
solid line corresponds to the Eddington limit. The dashed line separates our fiducial radio and quasar mode. The dotted lines corresponds 
to IG^^MEdtji 10~''^Eddi 10~*A^Edd) ^r"! 10~^MEdd from top to bottom. Upper left plot is for 2 = 0, upper right plot for z = 1, 
bottom loft plot for 2 = 2, and bottom right plot for 2 = 4. 



mode, and has more impact on the baryon content of galax- 
ies. 

We also tested the convergence and the robustness of 
our model to the effect of finite numerical resolution using 
different box sizes, DM particle masses, and minimum cell 
sizes. We obtain a satisfying convergence of the cosmic BIf 
density, and Mbu-Ms and i\fBH-o"s relationships at z = 
as long as the DM mass resolution is Mum ^ 4.4 10* Mq 
and the minimum physical cell size is Ax s; 1.52h"\kpc. 
However we have seen that the convergence of the cosmic 
SFR at high redshift is not reached because it is dominated 
by small halos. 

These simulations have demonstrated the ability of 
AGN feedback to efficiently suppress the amount of stars 
and cold gas. The removal of cold baryons begins at high 
redshift as soon as the first progenitors of massive galaxies 
collapse. At low redshift, the effect of AGN feedback is am- 
plified by the presence of more massive objects for which the 
impact of AGN feedback is stronger. The presence of AGN 



feedback is a necessary ingredient of galaxy mass budgets as 
it leads to a better fit to observed SFRs that slowly converge 
at low redshift. 

We have shown that quasars are increasingly important 
at high redshift. Since cold gas is more abundant in galaxies 
at high redshift, BHs accrete gas at higher rates, and, as 
a consequence, the quasar mode of AGN feedback is more 
often triggered. As the gas is consumed by star formation, 
and the accretion of new fresh gas is quenched by both the 
shock-heating of massive structures and the feedback from 
SN and AGN, the accretion onto BHs proceeds at lower 
rates. Thus, as time goes by, the radio mode of AGN feed- 
back becomes more and more dominant, and triggers jets in 
massive structures such as groups and clusters of galaxies. 
Radio active galaxies comprise only a fraction of observed 
AGN (Best et al. 2005; Smolcic et al. 2009, 2011). As the 
amount of radio emission varies a lot from object to object 
and depends on the gas accretion rate/mode onto the BH, 
working out the reasons of quantitative (dis)agreement be- 
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Figure 17. AGN luminosity functions at different redshifts for 
the 256L25JH simulation (left panels) and for the 256L50JH sim- 
ulation (right panels). The total luminosity functions (solid lines) 
are decomposed into the contribution from the radio mode (dot- 
ted lines) and from the quasar mode (dashed lines). 



tween our predictions and observations in detail, taking into 
account limitations on both sides, is a complex issue, beyond 
the scope of this paper. We therefore defer such an analysis 
to future work. 

Our results are in good agreement with previous cosmo- 
logical numerical simulations of the self-regulated growth 
of BHs. The quasar mode employed here is extremely 
similar to that modeled with a different numerical tech- 
nique by Booth & Schaye (2009), inspired by the model 
of Sijacki et al. (2007). Our conclusions in terms of Mbh /Ms 
relationships, and BH density evolution are, indeed, ex- 
tremely comparable. This is an important conclusion of this 
paper: cosmological codes that treat gravity and gas dy- 
namics with different approaches produce extremely similar 
results for the co-evolution of BHs and galaxies with an iden- 



tical set of sub-grid physics for galaxy formation. We have 
implemented a radio mode for AGN feedback and found a 
set of parameters that give results in good agreement with 
both observations and results produced by the well-tested 
quasar mode. It proves that a pure kinetic mode is also ca- 
pable of reducing the amount of baryons in galaxies and 
self-regulate the growth of BHs. A simple combination of 
both modes, the dual radio/quasar mode, reproduces the 
observations as well as any one of the two single modes, and 
provides a more realistic approach to the treatment of AGN 
feedback. 
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APPENDIX A: BULGE AND DISC 
DECOMPOSITION 

We decompose the bulge and the disc in our simulated galax- 
ies with a similar procedure to what is used for decomposing 
observed stellar luminosity profiles. The difference is that we 
do this decomposition directly on the stellar density on not 
on the luminosity. Our major assumption is that the rela- 
tion between stellar surface-luminosity and stellar surface- 
density is linear. The procedure is the following: 

- Galaxies with a minimum of 100 star particles are iden- 
tified with a galaxy finder based on the Most massive Sub- 
node Method described in Tweed et al. (2009), that allows 
for a clean separation of structures and sub-structures. This 
is particularly important for galaxy mergers, or for massive 
galaxies with extended stellar halos and galaxy satellites. 

- A stellar surface density profile is computed from the 
stars that belong to a given galaxy as defined by the galaxy 
finder, with the galaxy seen face-on. The face-on view of 
the galaxy is defined by the line-of-sight along the angular 
momentum axis defined by the rotation of the stars. 

- A double decreasing exponential of the form oc 
Eiexpr/r; (the "i" subscript is for the bulge Eb, rb, or 
for the disc Ed, rd) is fitted to the stellar surface-density 
profile using a least to separate the bulge from the disc 
component. 

- The bulge mass is taken to be Mh = 27rrbEb; the disc 
mass is taken to be Md = 27rrdEd. 

- If only one component fits the stellar density profile, we 
directly use the stellar mass given by the galaxy finder. 

We point out that, usually, cruder strategies are em- 
ployed to separate the so-called "bulge" component from 
a galaxy's total stellar mass, by assuming that the half- 
mass equals the bulge mass (Sijacki et al. 2007), or using 
the total halo stellar mass as a proxy for the bulge stellar 
mass (Di Matteo et al. 2008; Booth & Schaye 2009). Even 
though our method is more consistent with the observational 
definition of a bulge, our tests suggest that at these kpc res- 
olutions, other definitions lead to similar conclusions. 

Fig. Al shows the stellar surface-density profile ob- 
tained for the most massive galaxy in the 256L25JH simula- 
tion. We can clearly observe an inner bulge component with 
size rb = 8.12 kpc, and a disc component with larger charac- 
teristic radius rb = 21.85 kpc. This galaxy has a bulge mass 
Mb = 3.7110"Mq and a disc mass Md = 4.32 1O"M0. 
The sum of the two fits with the exponential forms is a good 
approximation of the stellar-density profile and gives a total 
stellar mass = 8.03 10^^ Mq similar to the total stellar 
mass obtained with the galaxy finder Mgai — 8.34 10^^ M©. 
We must stress that the disc component for such massive 
galaxies is not relevant, because this outer component of 
the galaxy corresponds to a large stellar halo rather than a 
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Figure Al. Stellar surface-density (black solid) for the most mas- 
sive galaxy at 2 = in the 256L25JH simulation with the bulge 
(red dotted) and disc (red dashed) decomposition from the best 
fitting model (red solid). 

rotating disc of stars. However, in this paper, we do not dis- 
cuss the properties of stellar discs or stellar halos of galaxies. 
Our main focus is to separate the bulge component from the 
total distribution of stars. 



